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1  Introduction 

Nucleate  pool  boiling  is  a  very  efficient  and  widely  used  cool¬ 
ing  method  with  a  wide  variety  of  heat  dissipation  applications  in 
industry  and  technology,  e.g.,  in  aircraft  thermal  management, 
electronics,  and  nanofluidics.  A  better  understanding  of  the  mech¬ 
anisms  and  details  of  the  multiscale  phenomena  involved  in  the 
process  of  boiling  is  needed  to  provide  an  insight  into  ways  to 
control  and  enhance  the  boiling  heat  transfer. 

Additives  have  long  been  recognized  and  studied  for  the  pur¬ 
poses  of  enhancing  the  boiling  heat  transfer  [1-4]  and  heat  trans¬ 
fer  with  surfactant  additives  in  pool  boiling  is  the  topic  of  active 
research  in  thermal  management  [5-10],  spray-cooling  [11], 
micro  and  nanofluidics  [12-16].  Additives  can  enhance  or  dimin¬ 
ish  the  effectiveness  of  boiling  heat  transfer  depending  on  their 
chemistry  or  concentration  [1,3,4,9,13-15,17-21],  On  a  macro¬ 
scopic  level,  solution  additives  contribute  to  dynamic  surface  ten¬ 
sion  [2,6,17,22]  and  modify  the  surface  wettability  [23].  Since 
nucleate  boiling  is  such  ubiquitous  thermal  management  method, 
there  is  a  sustained  interest  in  its  enhancing,  better  understanding, 
and  control  [5,11,12,16,24-27].  There  have  been  a  number  of 
microscopic  models  attempting  to  describe  nucleate  boiling  with 
additives  [28-30]  but  there  is  still  a  need  for  better  understanding 
on  a  molecular  level.  Recently,  the  effect  of  concentration,  chem¬ 
istry,  temperature,  and  pH  on  thermal  conductivity  and  viscosity 
of  bulk  aqueous  surfactant  solutions  were  studied  experimentally 
by  Zhou  et  al.  [27],  It  was  observed  that  while  the  addition  of  any 
surfactant  decreases  the  thermal  conductivity  of  the  solution,  non¬ 
ionic  surfactants  reduce  the  thermal  conductivity  more  than  the 
ionic  ones. 

Our  research  interest  is  in  the  advancement  of  molecular-level 
understanding  of  the  processes  involved  in  the  phenomenon  of 
boiling  with  additives.  The  present  effort  attempts  to  address  two 
specific  aspects  of  the  problem  of  interest:  exploring  the  effect  of 
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adding  the  surfactant  SDS  on  the  thermal  conductivity  of  the  solu¬ 
tions  at  room  and  boiling  temperature  as  well  as  the  concentration 
dependence  of  the  diffusion  of  SDS  in  water. 

We  carry  out  classical  molecular  dynamics  (MD)  simulations 
[31],  where  the  individual  atoms  are  approximated  by  spheres 
with  van  der  Waals  radii  and  partial  charges,  the  intramolecular 
interactions  are  approximated  via  (harmonic)  potentials  for  chemi¬ 
cal  bonds,  angles,  and  dihedrals,  while  the  intermolecular  interac¬ 
tions  are  modeled  by  the  long-range  Lennard-Jones  and  Coulomb 
interactions.  At  each  simulation  step,  the  Newtonian  equations  of 
motion  are  solved  i.e.,  numerically  integrated  for  each  atom  in  the 
potential  created  by  the  rest  of  the  atoms  in  the  system.  The  force 
field  parameters  are  usually  optimized  to  correctly  reproduce 
some  experimental  properties  and/or  are  developed  from  ab  initio 
or  density  functional  theory  (DFT)  calculations  of  gas  or 
condensed  phase. 

There  exist  a  large  number  of  theoretical  MD  model  potentials 
for  water  [32-38]  and  each  of  them  reproduces  well  only  some  of 
the  water  properties  and  only  within  certain  thermodynamic  con¬ 
ditions.  The  variety  of  water  models  can  be  broadly  divided  into: 
rigid,  flexible,  and  polarizable.  Rigid  models  keep  the  bond 
lengths  and  angles  at  their  equilibrium  positions  which 
significantly  reduces  the  computational  cost,  especially  for  large 
systems.  Additional  computational  savings  come  from  the  oppor¬ 
tunity  to  use  larger  simulation  steps.  Since  the  timestep  of  the 
integration  is  chosen  to  be  smaller  than  the  period  of  the  fastest 
motion  in  the  system,  usually,  a  timestep  of  1CT15  s  is  chosen  to 
account  for  the  fast  motion  of  the  hydrogen  atoms  in  water  mole¬ 
cules.  If  rigid  bond  models  are  used  via  algorithms  like  SHAKE 
[39]  the  timestep  can  be  increased  to  2  x  1CT15  s.  Flexible  water 
models  allow  the  water  molecule  to  deform,  i.e.,  bond  distances 
and  the  bond  angle  change  in  the  course  of  the  simulation,  at  the 
expense  of  increased  computational  cost.  Polarizable  force  fields 
additionally  allow  for  variable  partial  charges  on  the  atoms  of  the 
water  molecule  in  accordance  to  the  changing  environment  and 
this  is  why  they  are  exceptionally  computationally  costly.  All 
force  fields,  therefore,  account  for  the  many-body  effects  to  a 
varying  degree  since  including  the  polarization  effects  is  computa¬ 
tionally  demanding.  Moreover,  any  classic  force  field  description 
of  quantum  molecular  bonds  is  inherently  approximate  [40,41]. 
On  the  other  hand,  there  exist  a  variety  of  force  fields  that  are 
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parameterized  for  more  or  less  specific  system  or  class  of  systems 
in  mind.  For  example,  AMBER  [42],  CHARMM  [43],  and  OPLS 
[44]  force  fields,  to  mention  a  few,  are  optimized  for  biomolecular 
simulations,  while  CLAYFF  [45]  was  developed  with  the  goal  to 
be  general  enough  and  suitable  for  modeling  aqueous  solutions 
and  interfaces  of  inorganic  mineral  materials.  The  complexity  of 
the  problem  of  choosing  a  water  model  and  force  field  parameters 
to  describe  a  molecular  system  is  additionally  increased  by  the 
specific  suitability  and  compatibility  of  the  water  model  with  the 
model  parameters  to  describe  the  rest  of  the  simulated  system. 

Therefore,  the  choice  of  a  force  Held  depends  on  its  range  of 
applicability  and  on  the  molecular  system  to  be  studied  [34,46], 

Bulk  liquid  properties  using  rigid  water  models  are  extensively 
studied  and  available  in  the  literature  [32,34,38,41,47-49]. 

Recently,  the  transport  and  thermal  properties  of  well  established 
and  widely  used  rigid  and  flexible  water  models  were  studied  and 
compared.  Mao  and  Zhang  [46]  evaluated  the  thermal  conductiv¬ 
ity,  shear  viscosity,  and  specific  heat  of  seven  rigid  water  models. 

The  work  of  Sirk  et  al.  [50]  studied  the  thermal  conductivity  of 
common  classic  flexible  water  models,  compared  the  computa¬ 
tional  approaches  for  evaluating  thermal  conductivity,  and 
observed  that  although  both  rigid  and  flexible  water  models  over¬ 
estimate  the  thermal  conductivity  of  neat  water,  the  rigid  water 
models  provide  values  closer  to  experiment.  It  was  hypothesized 
that  the  increased  number  of  degrees  of  freedom  of  the  flexible 
water  models  is  responsible  for  the  overestimation  of  the  heat 
transfer.  Although  the  rigid  water  models  are  better  suited  for  sim¬ 
ulating  heat  transfer  in  bulk  fluids,  it  was  also  suggested  that  the 
appropriateness  of  rigid  versus  flexible  water  models  at  interfaces 
needs  further  investigation. 

We  report  a  computational  study  on  the  effect  of  surfactant  con¬ 
centration  on  the  diffusion  and  thermal  properties  of  aqueous  SDS 
solutions  at  room  and  at  boiling  temperatures.  We  also  assess  the 
performance  of  a  relatively  new  water  model  [38]  derived  from 
condensed  phase  ab  initio  calculations,  in  simulating  transport 
behavior  of  bulk  water  near  boiling  conditions  along  with  room 
temperature.  The  structure  of  the  paper  is  as  follows.  In  the  first 
part  of  the  paper,  a  model  for  SDS  by  Schweighofer  et  al.  [51] 
was  adopted  which  combines  all-atom  approach  for  the  hydro¬ 
philic  sulfate  headgroup  and  united-atom  course-grained  approach 
for  the  hydrophobic  tail  of  SDS.  The  diffusion  coefficient  of  SDS 
and  the  thermal  conductivity  of  the  aqueous  surfactant  solution 
are  computed  for  different  concentrations  of  SDS  at  room  and 
boiling  temperatures.  The  solvent  model  used  is  a  standard  TIP3P 
[49]  water  model.  The  results  for  the  diffusion  of  SDS  in  water 
and  the  thermal  conductivity  of  surfactant  aqueous  solutions  are 
compared  to  experimental  results.  In  the  second  part,  we  introduce 
the  ab  initio  flexible  water  model  developed  by  Akin-Ojo  et  al. 

[38]  in  2008  using  the  relatively  new  adaptive  force-matching 
method  [52,53]  and  compute  the  diffusion  and  thermal  conductiv¬ 
ity  of  bulk  water  predicted  by  this  model. 

2  Transport  Properties  of  Aqueous  Surfactant 
Solutions  at  Room  and  Boiling  Temperatures 

2.1  Surfactant  Model  and  Simulation  Details.  In  this  work, 
we  adopted  the  hybrid  model  of  SDS  [51,54—57]  which  combines 
all-atom  description  of  the  hydrophilic  sulfate  headgroup  of  SDS 
with  a  united-atom  approach  to  describe  the  hydrophobic  aliphatic 
tail  of  SDS,  bottom  structure  in  Fig.  1.  The  top  structure  of  Fig.  1 
shows  the  atomic  structure  of  the  anionic  surfactant  SDS.  An  all¬ 
atom  approach  describes  each  atom  with  its  own  set  of  parameters 
and  was  used  to  describe  the  water  molecules  in  this  work.  The 
united-atom  approach  is  in  a  sense  a  coarse-graining  approach  in 
which  a  group  of  atoms  is  described  as  one  "united-atom”  with  its 
own  set  of  parameters.  Only  the  hydrophobic  tail  is  described  by 
united-atom  approach  where  the  hydrogen  atoms  belonging  to 
each  carbon  atom  are  grouped  into  a  new  united-atom.  For 
example,  the  terminal  CH3  group  of  SDS  is  represented  by  the  C3 
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Fig.  1  Top-.  All-atom  model  of  SDS.  Bottom:  Hybrid  all-atom 
united-atom  model  of  SDS. 

atom  in  united-atom  approach.  The  CH2  groups  are  represented 
by  a  new  atom  type  CT,  and  the  last  CH2  group  bonded  to  the 
hydrophilic  head  is  designated  as  C2. 

2.1.1  Self-Diffusion  of  SDS  in  Water.  The  diffusion  coeffi¬ 
cient  of  a  solute  in  water  can  be  computed  from  the  mean  square 
displacement  (MSD)  of  the  center  of  mass  of  each  solute  molecule 
as  follows: 

D  =  ilim^{|r(/)-r(0)|2)  (1) 

6  r-too  at 

where  (|r(f)  —  r(0)|2}  is  the  mean  square  displacement  averaged 
over  all  solute  molecules  and  time  origins. 

In  the  simulations  for  obtaining  the  diffusion  coefficient  of  SDS 
in  water,  the  water  potential  used  was  TIP3P  [49],  Two  sets  of 
simulations  were  performed:  with  rigid  and  flexible  TIP3P,  at 
room  temperature,  298.15  K,  and  at  boiling  temperature,  373  K. 
The  initial  size  of  the  simulation  box  was  100  A  x  100  A  x  100  A 
containing  approx.  33,000  water  molecules.  The  initial  simulation 
boxes  were  built  using  Packmol  [58]  and  were  further  minimized 
and  equilibrated  at  the  respective  temperature.  Periodic  boundary 
conditions  were  imposed  in  all  three  directions  to  simulate  bulk 
properties  and  particle-particle  particle-mesh  method  as  imple¬ 
mented  in  Large-scale  Atomic/Molecular  Massively  Parallel  Sim¬ 
ulator  (LAMMPS)  [59,60],  also  related  to  Ewald  summation 
[31,61],  was  applied  to  correctly  account  for  the  long-range  elec¬ 
trostatic  interactions.  All  simulations  were  carried  out  using  the 
LAMMPS  [60],  The  input  files  for  LAMMPS  were  prepared  using 
VMD  [62].  After  a  standard  procedure  of  minimization  and  equili¬ 
bration,  we  performed  2-ns  equilibrium  MD  runs  at  constant  vol¬ 
ume  and  temperature  of  the  system  at  298.15  K  and  373  K, 
respectively,  in  order  to  obtain  the  self-diffusion  of  SDS  in  water. 

2.1.2  Thermal  Conductivity.  Thermal  conductivity  can  be 
computed  from  equilibrium  MD  simulations  using  the  Green- 
Kubo  formula  [63]  or  by  nonequilibrium  MD  (NEMD)  simula¬ 
tions  by  imposing  a  known  heat  flux  in  the  system.  We  computed 
the  thermal  conductivity  by  carrying  out  NEMD  simulations 
where  heat  flux  is  introduced  in  the  system  which  leads  to  estab¬ 
lishing  a  temperature  gradient  across  the  system  as  proposed  by 
Miiller-Plathe  [64],  We  used  the  Miiller-Plathe  method  of  impos¬ 
ing  heat  flux  [65]  which  involves  swapping  the  kinetic  energy  of 
the  hottest  particle  in  a  predefined  “cold”  region  of  the  system 
with  the  kinetic  energy  of  the  coldest  particle  in  the  predefined 
“hot”  region  of  the  system  thus  conserving  the  total  energy  of  the 
system  while  inducing  heat  flux  from  the  hot  region  (heat  source) 
to  the  cold  region  (heat  sink).  In  response  to  the  heat  flux  through 
the  system,  a  temperature  gradient  is  established.  After  equilibrat¬ 
ing  the  simulation  systems  at  the  respective  temperature  and  con¬ 
stant  volume,  NEMD  runs  at  constant  volume  and  energy  were 
performed,  during  which  a  steady  temperature  gradient  in  r-direc- 
tion  was  established.  From  the  last  700  ps  of  the  simulations,  the 
temperature  profiles  were  recovered  and  the  temperature  gradient, 
dT/dz,  was  computed  from  linear  fits  of  the  temperature  profiles. 
The  temperature  gradient  is  induced  by  the  heat  flux,  /,  imposed 


Table  1  Diffusion  coefficient  Ds  in  10-6  cm2/s  of  aqueous  SDS 
solutions  at  room  temperature,  298.15  K,  computed  for  a  range 
of  concentrations.  The  experimental  results  fl])  are  taken  from 
Ref.  [66], 


C  (mol/L) 

ExpU 

TIP3P(f) 

TIP3P(r) 

0.0016 

_ 

0.51 

8.18 

0.0050 

— 

7.16 

4.64 

0.0082 

— 

5.89 

4.81 

0.0100 

1.76 

2.87 

6.12 

0.0164 

3.00 

8.47 

12.34 

0.0500 

4.03 

6.32 

13.18 

0.1250 

4.53 

3.09 

4.23 

Table  2  Thermal  conductivity  l  (W/m  K)  of  aqueous  SDS  solu¬ 
tions  at  298  K  and  373  K  and  different  surfactant  concentration 


C  (mol/L) 

TIP3P(f)  TIP3P(r) 

298. 15K 

TIP3P(f) 

TIP3P(r) 

373  K 

0.0 

0.96 

0.82 

1.00 

0.68 

0.0016 

0.91 

0.80 

0.96 

0.83 

0.0082 

0.92 

0.72 

0.95 

0.76 

0.0164 

0.91 

0.80 

0.96 

0.84 

0.0500 

0.90 

0.80 

0.94 

0.83 

0.1250 

0.88 

0.79 

0.92 

0.82 

C  (mol/L) 


Fig.  2  Diffusion  coefficient  of  SDS  at  different  surfactant  con¬ 
centration  and  373  K  computed  with  flexible  TIP3P  (Flex)  and 
rigid  TIP3P  (Rigid)  water  potential 


on  the  system  and  the  coefficient  of  proportionality  is  the  thermal 
conductivity,  X,  of  the  material:  J  =  —XdT / dz . 

2.2  Diffusion  of  SDS  in  Water:  Results  and  Discussion. 

When  SDS  (NaC^^sStTp  dissolves  in  water,  the  ionic  bond 
between  Na+  and  DS~(Ci2H25SOj)  breaks  and  sodium  counter¬ 
ions  dissociate.  The  degree  of  dissociation  and  the  fraction  of 
sodium  ions  that  is  associated  in  the  first  and  second  solvation 
shells  can  be  determined  from  the  sodium-sulfur  radial  distribu¬ 
tion  function  [55].  For  the  purposes  of  computing  the  diffusion  of 
SDS  in  water  and  comparing  to  the  experimental  values  [66],  and 
since  the  light  sodium  ion  is  far  more  mobile  than  the  DS~  anion, 
hereafter  we  consider  the  diffusion  of  the  anion  DS~  and  refer  to  it 
as  SDS.  The  diffusion  coefficient  of  SDS  was  computed  from  the 
MSD  of  the  center  of  mass  of  each  molecule,  according  to  Eq.  (1). 
The  computed  values  at  298.15  K  are  shown  in  Table  1  and  com¬ 
pared  with  the  experimental  values  from  Ref.  [66].  It  should  be 
noted  that  the  computation  of  the  MSD  at  low  surfactant  concen¬ 
trations  is  inherently  less  accurate  than  at  larger  concentrations. 
Surfactant  concentration  of  0.0016  mol/L  corresponds  to  a  single 
SDS  molecule  in  our  simulation  box  whereas  a  surfactant  concen¬ 
tration  of  0.125  mol/L  corresponds  to  75  SDS  molecules  in  the 
simulation  box  which  provide  better  statistics  for  MSD.  From 
Table  1,  it  is  clear  that  the  computed  values  of  the  diffusion  coeffi¬ 
cient  of  SDS  in  water  are  of  the  correct  order  of  magnitude.  At 
higher  concentrations,  where  we  expect  better  accuracy,  the  rigid 
TIP3P  water  potentials  systematically  predicts  larger  values  than 
the  flexible  one.  Similar  trend  is  also  observed  at  boiling  tempera¬ 
ture,  as  it  can  be  seen  from  Fig.  2,  where  diffusion  coefficient  of 
SDS  computed  at  373  K  is  plotted  as  a  function  of  the  surfactant 


concentration.  Again,  the  diffusion  coefficients  of  SDS  in  rigid 
TIP3P  water  are  always  larger  than  the  values  predicted  by  using 
flexible  water  potential. 

The  calculations  from  simulations  with  rigid  water  model  seem 
to  systematically  overestimate  the  diffusion  coefficient  of  SDS  at 
higher  surfactant  concentrations.  At  lower  surfactant  concentra¬ 
tions,  the  accuracy  of  the  simulated  results  is  lower,  probably,  due 
to  insufficient  statistics  and  simulations  of  longer  duration  and/or 
larger  scale  can  improve  the  accuracy  of  the  estimates. 

2.3  Thermal  Conductivity  of  Aqueous  SDS  Solutions: 
Results  and  Discussion.  The  thermal  conductivities  of  aqueous 
SDS  solutions  at  a  range  of  surfactant  concentrations  were  com¬ 
puted  by  Miiller-Plathe’s  method  [64,65]  from  NEMD  simula¬ 
tions.  The  values  of  the  thermal  conductivity  of  aqueous  SDS 
solutions,  computed  with  flexible  and  rigid  TIP3P  water  model  at 
room  and  boiling  temperature  and  different  surfactant  concentra¬ 
tions  are  listed  in  Table  2. 

Recently,  an  experimental  investigation  of  the  effects  of  the 
surfactant  concentration,  temperature,  and  pH  on  the  thermal  con¬ 
ductivity  of  aqueous  solutions  of  nonionic,  cationic,  and  anionic 
surfactants,  was  reported  [27].  Compared  to  the  reported  values  in 
Ref.  [27],  both  water  models  lead  to  higher  than  experimental 
thermal  conductivity  of  the  aqueous  SDS  solutions.  As  in  the  case 
of  pure  water,  the  flexible  water  potential  leads  to  an  overestimate 
of  the  thermal  conductivity  of  the  solution.  Both  models  correctly 
predict  the  increase  of  the  thermal  conductivity  with  increase  of 
temperature.  Zhou  et  al.  [27]  experimentally  observed  an  increase 
of  the  room  temperature  thermal  conductivity  ratio  (X/XwMeI)  °f 
the  ionic  surfactants  SDS  and  SDBS  at  lower  concentrations, 
approx.  0. 1  wt.  %  for  SDS,  and  gradual  decrease  of  the  thermal 
conductivity  ratio  of  the  solution  with  increasing  the  concentra¬ 
tion.  At  lower  concentrations,  it  was  observed  that  the  thermal 
conductivity  of  the  ionic  surfactant  solution  was  larger  than  the 
thermal  conductivity  of  neat  water.  Our  simulations  did  not  repro¬ 
duce  the  enhancement  of  the  thermal  conductivity  of  the  surfac¬ 
tant  solutions  at  low  concentrations  but  did  reproduce  the 
experimentally  observed  slow  decrease  trend  at  concentrations 
higher  than  the  critical  micelle  concentration  (CMC).  Clearly,  fur¬ 
ther  simulations  on  a  much  larger  scale  are  needed  to  better  under¬ 
stand  the  thermal  properties  of  low  concentration  surfactant 
solutions,  as  well  as  to  explore  the  effects  of  the  chosen  surfactant 
and  solvent  models,  simulation  length  and  system  size  on  the 
computational  results.  Studies  with  increasingly  large  simulation 
systems  can  be  the  subject  of  future  works. 

3  MP2f  Water  Model  and  Simulation  Details 

With  respect  to  reproducing  the  transport  properties  of  water, 
and  more  specifically,  the  self-diffusion  of  water,  the  modified 
TIP3P  [49]  water  potential  used  in  this  work  is  the  most  com¬ 
monly  used  and  reliable  force  field  as  it  was  optimized  to  better 
reproduce  not  only  the  experimental  water  density,  heat  of 
vaporization,  and  dielectric  constant  but  also  the  experimental 
bulk  diffusion  coefficient  of  water. 
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The  flexible  force  field  MP2f  [38]  was  developed  by  the  rela¬ 
tively  new  adaptive  force  matching  method  from  condensed  phase 
quantum-mechanics/molecular-mechanics  (QM/MM)  calcula¬ 
tions.  The  force  matching  (FM)  method  was  first  proposed  in 
1994  [52].  The  FM  method  optimizes  the  force  field  parameters 
not  to  best  reproduce  any  experimental  bulk  water  properties,  as  is 
the  case  in  the  empirical  force  fields,  nor  to  best  fit  single  point 
energies  from  ab  initio  gas  phase  calculations,  as  is  the  case  in  the 
usual  ab  initio  force  fields,  but,  instead,  to  reproduce  ab  initio 
forces  derived  from  DFT  calculations  in  the  condensed  phase 
[67].  The  adaptive  force  matching  (AFM)  method  is  an  improve¬ 
ment  upon  the  original  FM  schemes  in  that  it  uses  QM/MM  calcu¬ 
lations  and  avoids  the  need  for  expensive  condensed  phase 
electronic  structure  calculations  [38].  Therefore,  one  major  appeal 
of  the  AFM  method  lies  in  the  fact  that  the  force  field  parameters 
are  derived  from  condensed  phase  calculations  as  opposed  to  the 
usual  gas  phase  calculations,  since  those  force  field  parameters  are 
intended  to  be  used  to  simulate  bulk  liquid  water. 

In  order  to  investigate  the  transport  properties  of  the  MP2f 
water  model,  we  computed  the  diffusion  constant  of  the  model 
and  the  thermal  conductivity  of  water  at  room  temperature 
298. 15  K  and  boiling  temperature  373  K. 


Fig.  3  Diffusion  coefficient  as  a  function  of  the  inverse  length 
of  the  periodic  simulation  box.  From  a  linear  fit  of  the  data,  the 
size-independent  diffusion  coefficients  for  MP2f  water  model 
were  extrapolated  to  be  4.06  xIO-5  cm2/s  at  298.15  K  and 
9.95  x  10-5  cm2/s  at  373  K. 


3.1  Self-Diffusion  of  Bulk  Water.  We  performed  2-ns  equi¬ 
librium  MD  runs  to  compute  the  diffusivity  of  water  at  two  tem¬ 
peratures,  298.15  K  and  373  K,  using  MP2f  water  potential.  The 
simulation  step  was  0.001  ps.  The  simulation  runs  were  performed 
at  constant  volume  and  temperature  using  Nose-Hoover  thermo¬ 
stat  as  implemented  in  LAMMPS  [60].  Snapshots  of  the  atomic 
positions  and  velocities  were  output  every  100  simulation  steps  or 
every  0.1  ps.  The  self-diffusion  coefficient  of  water  is  computed 
from  the  MSD  of  the  center  of  mass  of  each  water  molecule  or,  to 
a  very  good  approximation,  the  position  of  the  oxygen  atom  in 
each  water  molecule,  according  to  Eq.  (1). 

As  it  is  well  known  that  the  computed  values  of  the  diffusion 
coefficient  depend  on  the  size  of  the  simulation  box  [68-70]  when 
periodic  boundary  conditions  are  applied,  we  set  up  a  series  of  six 
simulation  systems,  containing  128,  256,  512,  1024,  2048,  and 
4096  water  molecules,  at  normal  water  density  of  33.0  nm-3.  If  L 
is  the  length  of  the  periodic  water  box,  and  DL  is  the  computed 
diffusion  coefficient  in  the  same  periodic  water  box,  then  the  de¬ 
pendence  of  Dl  on  the  length  L  is  given  by 


Table  3  Size-independent  diffusion  coefficient  ( D0 )  and  viscos¬ 
ity  (//)  of  MP2f  water  model  at  298.15  K  and  373  K.  For  compari¬ 
son,  estimates  of  the  size-independent  diffusivity  of  other  water 
models  and  experimental  values  are  also  provided.  (*)  Ref.  [73], 
(*)  Ref.  [74],  (s)  Ref.  [69],  (*)  Ref.  [70],  (b)  Ref.  [46],  (#)  Ref.  [46], 
this  value  was  calculated  at  363  K. 


Model 

D0[x  10" 

"5cm2/s] 

ij[x  10 

4  kg  (m  ■  s)] 

298.15  K 

373  K 

298.15  K 

373  K 

MP2f 

4.06 

9.95 

3.24 

3.84 

IP3P 

6.05§ 

— 

3.18b 

2.17# 

SPC/E 

2.97* 

— 

6.4* 

— 

TIP4P/2005 

2.49* 

— 

8.3* 

— 

Exp. 

2.3* 

— 

8.96# 

— 

Table  4  Computed  thermal  conductivity  l  in  W/(m  K)  of  neat 
water  at  room  and  boiling  temperature.  The  experimental  values 
are:  (s)  from  Ref.  [75]  and  (+)  from  Ref.  [76]. 


Dl  =  Do  — 


2.831kBT 
6  ni]L 


(2) 


In  Eq.  (2),  r]  is  the  shear  viscosity  of  water.  Therefore,  we  can 
extrapolate  the  size-independent  diffusion  coefficient  D0  and  the 
viscosity  of  water  from  the  series  of  simulations  of  increasing  sim¬ 
ulation  box  size. 


298.15  K 


373  K 


TIP3P(f) 

0.96 

1.00 

MP2f 

0.64 

0.66 

TIP3P(r) 

0.82 

0.68 

Exp. 

0.607s 

0.6723 

3.2  Self-Diffusion  and  Thermal  Conductivity  of  MP2f 
Water  Model:  Results  and  Discussion.  The  values  of  the  self¬ 
diffusion  coefficient  of  water  obtained  from  simulations  in 
different  size  periodic  simulation  boxes  were  fitted  with  a  linear 
function  of  \/L  as  shown  in  Fig.  3.  From  these  fits,  the  size- 
independent  values  of  the  diffusion  coefficients  in  infinitely  large 
bulk  system  were  extrapolated.  The  fitting  procedure  provides  an 
estimate  of  the  viscosity  of  the  water  model  as  well.  The  extrapo¬ 
lated  size-independent  values  of  the  diffusion  coefficient  of  MP2f 
water  model  are  shown  in  Table  3.  In  the  table,  we  provide  also 
literature  values  obtained  by  the  same  extrapolation  technique. 
The  three  literature  water  models  are  rigid.  Our  estimate  of  the 
viscosity  of  the  MP2f  water  models  is:  3.24  x  lO^kg/fm-s)  at 
298.15  K  and  3.84  x  10^kg/(m-s)  at  373  K.  From  Table  3,  it  can 
be  observed  that  MP2f  water  potential  is  somewhat  better  than 
rigid  TIP3P  with  respect  to  the  predicted  value  of  the  size- 
independent  diffusion  coefficient  but  SPC/E  [71]  and  TIP4P/2005 
[72]  predict  size-independent  values  of  diffusion  coefficient  and 


viscosity  of  water  in  better  agreement  with  experiment.  Table  3 
also  lists  recently  reported  values  for  the  shear  viscosities  of  a 
number  of  rigid  water  models  computed  using  the  Green-Kubo 
method  [46]. 

The  values  of  the  computed  thermal  conductivity  of  water  are 
shown  in  Table  4.  It  was  reported  [50]  that  the  computed  value  of 
the  water  thermal  conductivity  was  found  to  be  independent  on 
the  size  of  the  simulation  box  (provided  it  is  larger  than  40  A), 
therefore,  we  performed  all  simulations  in  a  periodic  box  contain¬ 
ing  4096  water  molecules.  The  results  in  Table  4  confirm  the 
observation  from  Ref.  [50]  that  the  flexible  water  models  tend  to 
overestimate  the  thermal  conductivity  of  water,  TIP3P(f)  versus 
TIP3P(r).  Our  results  further  show  that  the  flexible  MP2f  water 
model  performs  exceptionally  well  in  reproducing  thermal  con¬ 
ductivity  of  water  at  room  and  boiling  temperature,  compared  to 
both  the  rigid  and  flexible  TIP3P. 

In  conclusion,  the  flexible  MP2f  water  potential,  derived  by 
AFM  method  from  condensed  phase  QM/MM  calculations  is 
shown  to  perform  reasonably  well  in  reproducing  the  diffusion 
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coefficient  of  water  and  very  well  in  reproducing  the  thermal  con¬ 
ductivity  of  water  at  both  room  and  boiling  temperature. 

4  Conclusions 

We  performed  an  extensive  series  of  MD  simulations  comput¬ 
ing  the  transport  properties  of  bulk  aqueous  surfactant  solutions  at 
room  and  boiling  temperatures  for  a  range  of  surfactant  concentra¬ 
tions.  The  effect  of  the  water  model  on  reproducing  the  surfactant 
diffusivity  and  the  thermal  conductivity  of  aqueous  SDS  solutions 
was  investigated  by  using  flexible  and  rigid  TIP3P  water.  It  was 
observed  that  at  higher  surfactant  concentrations,  the  rigid  water 
model  overestimates  the  values  for  surfactant  diffusion.  At  lower 
concentrations,  the  results  were  less  reliable  due  to  insufficient 
statistics  and  further  simulations  with  larger  systems  are  required 
to  improve  the  accuracy  of  the  estimates.  It  was  recently  observed 
in  the  case  of  pure  water  [50],  that  the  flexible  water  models  are 
not  only  more  computationally  costly  but  also  they  tend  to  overes¬ 
timate  the  thermal  conductivity  of  neat  water.  Our  results  for  the 
thermal  conductivity  of  aqueous  SDS  solutions  confirm  this  but 
also  show  that  when  the  thermal  conductivity  enhancement  of  the 
solution  with  respect  to  pure  water  is  concerned,  the  flexible  water 
model  reproduces  the  experimental  trend  more  reliably,  at  least  at 
concentrations  higher  than  CMC.  Future  larger-scale  MD  studies 
are  needed  to  further  explore  and  better  understand  the  transport 
properties  of  surfactant  solutions  at  low  concentrations. 

We  also  studied  the  ability  of  MP2f,  an  ab  initio  water  potential 
derived  from  condensed  phase  calculations,  to  reproduce  the  dif¬ 
fusion  coefficient  and  thermal  conductivity  of  water  at  room  and 
boiling  temperatures.  The  size-independent  diffusion  coefficient 
of  water  was  derived  from  a  series  of  MD  simulations  with 
increasing  size  of  the  simulation  system.  The  bulk  water  self¬ 
diffusion  coefficient  was  estimated  to  be  4.05  x  10~5cm2/s  at 
298  K  and  9.95  x  10~5cm2/s  at  373  K.  The  computed  values  of  the 
thermal  conductivity  were:  0.64W/(m  K)  at  298  K  and  0.66  W/ 
(m-K)  at  373  K,  in  much  better  agreement  with  the  experimental 
values  compared  to  both  the  rigid  and  the  flexible  TIP3P. 
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